An insight into the electronic structure of graphene: from monolayer to multi-layer 
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i In this paper, we analytically investigate the electronic structure of Bernal stacking (AB stack- 

ing) graphene evolving from monolayer (a zero-gap semiconductor with a linear Dirac-like spectrum 
around the Fermi energy) to multi-layer (semi-metal bulk graphite). We firstly derive a real space 
analytical expression for the free Green's function (propagator) of multi-layer graphene based on 
the effective-mass approximation. The simulation results exhibit highly spatial anisotropy with 
three- fold rotational symmetry. By combining with the STM measurement of d 2 I/dV 2 (the sec- 
ond derivative of current), we also provide a clear high-throughput and non-destructive method to 
identify graphene layers. Such a method is lacking in the emerging graphene research. 

PACS numbers: 73.61Wp, 61.72.Ji, 68.37.Ef 
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rjrj ' I. INTRODUCTION 

<D . 

The successful synthesis of monolayer graphite (graphene) has attracted enormous research interest over the past 
years.— Researchers can also experimentally manipulate few layer within multi-layer graphene samples and observe 
their quasi two-dimensional behavior £2d£&2. The conduction band of graphene is well described by the tight-binding 
model, which includes the it orbitals that is perpendicular to the plane each carbon atom locates. This model describes 
1 a semi-metal, with zero density of states at the Fermi energy. Here, the Fermi surface is reduced to two inequivalent 
K-points located at the corners of the hexagonal Brillouin Zone. The low-energy excitations with momenta in the 
q ■ vicinity of any Fermi points have a linear dispersion. They can be described by a continuous model, which reduces to 
O \ the Dirac equation in two dimensions. These observations have been validated experimentally 1 . 

Recent research attention has turned to study multi-layer graphene^^&i, because its electronic proprieties are 
I ' quite different from those of the monolayer graphene. Two-dimensional (2D) multi-layer graphene is an intermediate 
, crystals between bulk graphite and a single graphene plane. Relatively weak inter-layer coupling within multilayer 
graphene inherits some properties from monolayer graphenes. The special geometry of the Bernal stacking (A-B stack- 
ing) between graphene layers results in low-energy states mainly reside around B site. In contrast to corresponding 
three-dimensional (3D) bulk structures, electrons in a multi-layer graphene are confined along one crystallographic di- 
rection, which offers its unique electronic characteristics. From the Hall effect measurements of a multi-layer graphene, 
the results show that multi-layer graphene behaves like multi-carriers (coexistence of electrons and holes) semi-metallic 
systems*^ Moreover, carriers can switch from electron to hole transport by alteranting gate voltage. This phenomenon 
makes multi-layer graphene a remarkable platform for studying mesoscopic transport in low-dimensional materials. 
It is also worth noting that multi-layer graphene systems may be of interest for buidling nanoscale devices because 
we can easily change their electronic properties by using conventional lithographic techniques. 

To date, few studies attempting to uncover the unique electronic properties of multi-layer graphene in an analytical 
fashion. Moreover, the physical properties of multi-layer graphene are closely related to its propagator. In this 
paper, we first develop an analytical formula of the free Green's function for electrons in multi-layer graphene (Bernal 
stacking) in real space by using the effective-mass approximation. We then calculate the LDOS of multi-layer graphene 
and isotropic function by using Green's function explicitly with increasing number of graphene layers. Our numerical 
results show that monolayer graphene is dramatically different from other multi-layer graphene, and bilayer graphene 
capture the main characters of the multi-layer graphene. Currently only AFM and Raman spectrum ^ 10 ' 11 are used 
to detect the layers of graphene. In this paper, we present an effective method to detect graphene layers based on the 
STM. According to our derivative of LDOS with energy (proportional to the second derivative of current cPl/dV 2 ), 
we can quickly, clearly and nondistructively identify the graphene layers of finite multi-layer graphene based on STM 
images. 
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II. TIGHT-BINDING DESCRIPTION 



Figure 1 shows schematically the crystal structure of four AB-stacked graphene layers and the corresponding first 
Brillouin zone with the labels of the symmetric points. It was shown in the reference [12] that the interplanar spacing 



2 




PL 



B 



FIG. 1: (color online) The crystal structure of four AB stacked graphene layers. Type- A atoms (gray) have direct neighbors 
in the adjacent atomic layers while type-B atoms (yellow) are above hollow sites, (a) Side view with the corresponding tight- 
binding parameters, (b) Top view, (c) The first Brillouin zone of a finite number of graphene layers with the labels of special 
symmetric points (Here, we choose four-layer graphene as an example to illustrate our design). 



parameter a and lattice spacing d in two-layer graphene are almost identical to that in bulk graphite. We, therefore, 
choose the bulk value a = 1.42A and d — 3.35A in this paper. In order to study how the electronic structure 
of graphene changes from monolayer into multi-layer, we construct the tight-binding Hamiltonian for an arbitrary 
number of graphene layers. We extend our previous studies of monolayer graphene^ and of bilayer graphene^ to 
multi-layer one. Starting from two-layer graphene, type A and B carbon atoms are inequivalent: type-A carbon atom 
has a direct neighbor in its adjacent layer while type-B carbon atom does not and locates on the hollow sites. Since 
the physical properties of graphene are determined by tt bands near the Dirac point, only the contribution from the ir 
band is considered in this paper. In our tight-binding model, we limit our calculation to the nearest interaction carbon 
atoms on the graphene's p z orbitals. These interactions include interactions between the nearest A-B carbon atoms 
within a plane, interaction between nearest A-A carbon atoms of two nearest-neighbor planes. The tight-binding 
Hamiltonian of the multi-layer graphene is: 



where (r\i l ) is the wavefunction at site i on layer I. L is the number of graphene layers. V\ is the nearest hopping 
parameter within the layer and V2 is the nearest hopping parameter between two layers. We select V\ — — S.OeV and 
V% = OAeV in our calculation. 

The Bloch orbits for two nonequivalent sites, A and B, on layer I are written as |k^) = — = ^ e lk ' r 'A \ri A ) , 



|k' B ) = ^= J>2 e lk ' t 'b \ri B ). The summation is over all sites A and B on layer I. ri A and r/ B denote the real-space 

Ib 

coordinates of A and B sites, respectively. Here, N is the number of unit cells in the crystal. The tight-binding 
Hamiltonian can be obtained from the following matrix elements: 



L 




(1) 
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(k l A \H\k l A ) = (k l B \H\k l B )=0, 
/ Vifj,* (I = odd) 



<y A \H\k l B ) - , ., 
N A ' ' al \ Vi/j, (I = even) 

(k l B \H\k A ) = [V* (l = odd) 
x Bl 1 Al \ Vifj, (l = even) 

(k^lk^ 1 ) = {kT\H\k l A )^V 2 . 



(2) 



with 



y/3ak x ak y \f%ak x sL,, 
H = exp(ik y a) + exp[t( — )] + exp[i{ — — )]. 

In these caculations, we neglect all overlap integrals because of the relatively large separation between the carbon 
atoms. Furthermore, we are only interested in small energy electronic structure around the K point. With the 
expressions in Eq (2) , we can construct the tight-binding Hamiltonian for an arbitrary number of layers 



/ V lU * V 2 • 

V\u • 

H= v 2 o o vm ■ 

Vifj,* • 

V ; ; ; : 

The eigenfunction for multi-layer graphene system is now given by 

L L 



(3) 



|k)=^^|k5 4 )+^^ B |kL). 



(4) 



i=i i=i 

The hopping between sublattices of different layer lead to the coupled Harper equations, 



= v^x 1 + Vin*<t> B + v^X 1 {l = odd) 

£ (j) l B = Viu4> l A 



and 



e4> l A = V 2 ^ A L +V 1 ^ 
e<p l B = V 1 n*4 



A (I = even). 



(5) 



(6) 



where I = 1, 2, • • • , L denotes the layer number. The open boundary condition in the direction perpendicular to the 2D 
graphene plane is <fP A — 4> A +1 — 4>% = 4> B +1 = 0. According to the coupled Harper equations Eq.(5) and Eq.(6), the 
wave functions of the eigenstates and the corresponding energy spectrum can be obtained in the following analytical 
form. 



(7) 



b A (m, s,k) = D(m,s,k)sin( 
b l B (m,s,k) = 



mlir 
L + l 

Viti 



), 



D(m, S ,k)^sm(f£) (l = odd) 
D(m,s,k)^ Sm (^) (I = even) 



(8) 



where D(m, s, k) = { 53 [1 • 
i=i 



^^}sin 2 {^)}- 1 / 2 is the normalized factor (to = 1, 2, 



, L) and s = ±1. 



4 



III. REAL-SPACE GREEN'S FUNCTION 

Based on the definition of retarded Green's function, we obtain the reciprocal space Green's function for multi-layer 
graphene as: 



G$ ,,(m,k,£) = 



s=±l 



^("t,5,k)^'(m,a,k) 



(9) 



where i] is an infinitesimal quantity, /i and v label A or B. By taking the Fourier transform of G r e l^, (m, k, E) in the 
first Brillouin zone (1BZ), we can obtain the exact expression of real space multi-layer graphene Green's function, 



G-'dvr' iE ) = 



L 

E 

rn—l 



dkG r °*(m,k,E)e 



ik-(r I(J -rJ, ) 



(10) 
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Based on the effective-mass approximation, the integral is carried out around six corners in the first Brillouin zone 
within the low-energy region, which can form two 360° integrals around two incquivalcnt Dirac points. This method 
has been well used to study the electronic structure of monolayer and bilayer graphene in our previous works i 13 ' 14 
The real space multi-layer Green's function can be obtained as following: 



G r et (r lA ,v 
G r o e \r lB ,v 

Gl et {v lB ,t 



,,E) = cos\K 1 -(r u -r' l ,)]F 1 (\r u -r' ll \,E), 



coslK 1 ■ (ri B — v' v )]F2(|ri B — rj, \,E) (I — odd, I' — odd or I — even, I' — even) 



,, E ) = { 

,E) = 
,E) = 



cosiK 1 ■ {r lB - v' v ) - 2a r y, ]F 3 (\ri B - v[, \,E) (I = odd, I' 



even) 



where 



Fi(\r, 



F 2 (K 



F 3 (\ri, 



1 1' I ) 



E) 



• E) 



W' I) 



H' I) 



, }F 3 (\r lB - r' v \,E) (I = even,!' = odd) 

B B 

}F 4 (\r lA -r' l , \,E) (V = odd) 



cosfK 1 ■ (ti b — t' v ) + 2a r 
sin\K l ■ (r tA - r',, ) + a r , 

B " l B 

sin[K ■ (r tA - r' v ) - a r y ]F 4 (|r u - v\, \,E) {V = even) ' 

B 'b b 

smpK 1 ■ (ri B - r' v ) - a r y ]F 4 (|r iB - r' v \, E) (I = odd) 

A 

sin[K ■ {r lB -r' v ) + a r y }F 4 (\r lB - r' \,E) {I = even) ' 



_ kb\ m ,s,\t) S in{^) S in(^)J {k\v^-r\,J\) 



^ ^ — / ,//,• 



EE — 

s— ± 1 m— 1 



7T 



dk 



k 3 D 2 ( m ,sM)sin(fj^)sin(^i T f)Jo(k\r l ^-r' l ,J) 

{ V 2Cos( ^) +s[V 2 cos 2 { J^_ ) + {rk) 2 ] l/2 } 2 

E + in- {V 2 cos(^) + s[V 2 W(f^) + (rt) 2 ] 1/2 ) ' 



E) = 



E) 



L 

EE 

s— ± 1 m— 1 



2S 7 2 



dk 



EE — 



dk 



fc 3 D 2 (m,s,k)sm(fii^)sm(^ F f )J 2 (fc|r, M -r;, |) 

{V 2 COs(f^) + s[V2 C Os2(f^) + (rk)2Y/2 } 2 

E + in- OWf^) + s [V 2 cos 2 (f^) + (rfc)2]Va> ' 

k 2 D 2 ( m , S M)sin(f^) S in(^+).J 1 (k\r l ^-r' l ,J) 
V 2 cos(^)+s[V 2 co S 2 (^f T ) + (rk) 2 ]^ 2 " 



t^f^ * Jo E + i v -{V 2 cos(^) + 8^008*^) + (rk)*]V*Y 



(11) 



(12) 



D(m, s,k) is the normalized factor under the effective mass approximation. It can be obtained from the expression 
of D(m, s,k) by substituting V 2 \x[i* with (jk) 2 . Here, 7 = ZaV\/2. a rii y u is the angle between r M — r' u and x axis. 

K 1 = (47r/3\/3a, 0) and J n is the type-I n-order Bessel function. S — 3\/3a 2 /2 is the area of the unit cell in real 
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FIG. 2: The band structure of (a) monolayer, (b) bilayer, (c) trilayer and (d) 4-layer bernal stacking graphene along YKMY. 
The inset is an enlargement of the region indicated by the rectangular around the K point. 
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space and k c is the cutoff wave vector. From Eq.(ll) and Eq.(12), it is clear that the real space Green's function 
of multi-layer graphcne is constructed by multiplying two terms: spatially anisotropic and spatially isotropic terms. 
The spatially anisotropic (three-fold rotational symmetry) can be represented by sine and cosine functions which can 
be determined from two nonequivalent sets of the Dirac points. The second term including F%, F 2l F3, F4 is spatially 
isotropic in real space and depends only on the distance between two sites. For L=l and L=2, the expression of the 
Green's function obtain in Eq.(ll) and Eq.(12) can directly reduce to our previous results for monolayer and bilayer 
graphenei 13 i 14 

IV. NUMERICAL RESULTS 
A. Band structure 

As an example of applying the dispersion relation obtained in Eq.(7), we compare the energy band structure of 
bernal stacking multi-layer graphene involves from monolayer to 4-layer, the corresponding energy bands along the 
TKMT lines in the first Brillouin zone are shown in Fig. 2(a) ~ 2(d), respectively. In all the calculations, the Fermi 
energy is set to Ej = 0. The following are our observations: 

(i) For a monolayer graphene structure as shown in Fig. 2(a), the two bands cross at the K point leads to the fact 
that a monolayer graphene is a zero-gap metal. Around the K point the spectrum is linear as shown in the inset and 
it is given by ef = sjk. Here, k is wave vector from the Dirac K point. Comparing this expression with the relativistic 
energy expression e = \J m 2 c^ + p 2 c 2 , one can see that the dispersion relation of monolayer graphene mimics a system 
of relativistic Dirac particles with zero rest mass and an effective speed of light c = 10 6 m/s, which is 300 times smaller 
than the speed of light in vacuum. 

(ii) For a bilayer graphene structure as shown in Fig. 2(b), the number of levels is doubled. The spectrum is clearly no 
longer linear around the K point, but is parabolic, ef = V 2 /2 + sy/(V 2 /2) 2 + (■yk) 2 , ef = -V 2 /2 + s v / (V 2 /2) 2 + (yk) 2 . 
If moving away from the K point, the spectrum becomes linear again. In the present calculation by only considering 
the nearest interaction between layers, we observe that the bands just touch at the K point. However, a more detailed 
investigation based on first principle method shows a small overlap and an interaction leading to anticrossings between 
the conduction and valence bands^ 

(iii) For a trilayer graphene structure as shown in Fig. 2(c), there are four bands cross at the K point around zero en- 
ergy, two of them are linear, ef = V2V 2 / 2 + s {y/2V 2 / 2) 2 + ( 7 fc) 2 , ef = sjk, ef = -V2V 2 /2 + s^{V2V 2 /2) 2 + (7/c) 2 . 
It is clear that the band structure around the K point becomes more complex with increasing number of graphene 
layers than monolyaer and bilayer graphenes. For instance, the number of layers around the Fermi energy at the K 
point is doubled in comparison with monolayer and bilayer graphenes. The band diagram can be understood as the 
combination of band structure in monolayer graphene and in bilayer graphene. 

(iv) For a 4-layer graphene, shown in Fig. 2(d), four subbands cross exist at the K point. From Eq.(7), in the frame 
work of the nearest tight-binding model a general number of the subband cross at the K point can be deduced. For 
the number of layers (L), there are 2L subband. L of them cross at the K point if L is even, while L+l subbands cross 
at the K point if L is odd. The additional one subband, existing only when cos(j^) = (i.e. m = ^5-^), is linear. 

B. Surface Green's Function 

In this section, we study the behavior of the top layer Green's function with the changing number of graphene layers 
(L) in detail. From Eq.(ll) and Eq.(12), the corresponding real space Green's function of the top layer (I = I' = 1) 
are 

co8\K 1 -(ti a -t> 1a )]Fi(\t 1a -i' 1 J,E) 
COS [KMri B -riJ]F 2 (|r lB -r' lB |,£) 
siniK 1 • (r u - r' lfi ) + a^^jF^r^ - r[J,E) 
sinlK 1 ■ (r lB - r' u ) - a riB , r ^]F 4 (|r lj3 - r[ A \,E) 

(13) 

From previous studies, we find that in the low bias voltage, sites A and B are equivalent for monolayer graphene. 
They form honeycomb lattice in the STM image. While for a bilayer graphene, sites A and B are inequivalent, and 



G r et (r lA ,r' lA ,E) = 

G^{v lB ,v' lB ,E) = 

G r et (r lA ,r' lB ,E) = 

G r *(r lB ,r[ A ,E) = 
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FIG. 3: (color online) (a) and (b) are the LDOS of A and B site at the top layer with different graphene layers (from down to 
up L=l ~ 7), (c) and (d) are the corresponding derivative of LDOS with energy. To clearly see the difference between them, 
we shift each curve in (a) ~ (d) upwards. 



the LDOS is residual at site B and they form triangular lattice. Therefore, one direct approach to find the difference 
among multi-layer graphene is to study the LDOS at sites A and B by changing the number of stacked layers. From 
Eq.(I3), LDOS at sites A and B at the top layer can be expressed as: 
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FIG. 4: (color online) Spatially isotropic function F\, F2 and F4 of the top layer with different graphene layers (L=l ~ 7). 
Here we set E = O.leV. 
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Figure 3(a) and (b) present our numerical results of LDOS at sites A and B with different number of graphene 
layers. The LDOS for L=l is linear, which is clearly different from other curves. Starting from h—2, increasing the 
number of stacked layer, the LDOS just shows a slight deviation from its bilayer value for both A and B sites. This 
character means that monolayer graphene is a special case in the family of multi-layer graphene. The relatively weak 
inter-layer coupling within multi-layer graphene not only inherits some properties from monolayer graphenes, but also 
solidify some proprieties between different stacked layers. Detailed studies also find that LDOS at A site converge 
more slowly than B site with increasing the graphene layers (Twenty layers are needed for the LDOS at A site to 
converge to its bulk value, while just ten layers needed for B site). LDOS at A site fluctuate more dramatically than 
that at B site in the process of convergence. 

Since the first synthesis of monolayer graphene, except AFM and Raman spectrum ^ 10 ! 11 , few method can be used 
to exactly detect the number of graphene layers. In our work, based on STM second derivative map, we develop a new 
method to measure graphene layers, which is lacking in the emerging graphene research. From Eq. (14), the derivation 
of the LDOS with respect to the energy (dpo/dE) is proportional to the second derivation of current (cPl/dV 2 ) in 
the STM measurement, which corresponds to the excitation of vibrational modes 16 . The results are shown in Fig. 3(c) 
and (d). The curves are symmetry to the positive and negative energy. So for simplicity, we just consider the positive 
energy section in the following discussion. For site A, monolayer graphene has no peak, bilayer has one and trilayer 
has one as well. General relation for the peak number equals Int(L/2) (the function Int round off the variable to 
an integer). For L is even number, L and L+l has the same number of peaks, but the peaks for L+l move towards 
the high energy direction relative to the L. With increasing the number of stacked layers, the intensity of the peak 
become much smaller. For L—2, the peak height is about 0.16 eV~ 2 . While for L=14, the highest peak reduce to 0.06 
eV~ 2 . Additionally, the position of small peak is hard to determine when L become large than 14. Therefore, in the 
real experiment, only the sample has a few layers (less than 14 layers) may be detected by considering the sensitivity 
of the STM measurement. For site B, we observe no peak exists, but just step curve. With increasing stack layer, the 
step become not so clearly anymore. 

In the above discussion, we mainly focus on the LDOS of the top layer for the multi-layer graphene. However, 
to study the propagating behavior of electrons in the top layer of multi-layer graphene, we also need to know the 
properties of the space isotropic function F\, F 2 , and F4. in Eq.(12). The numerical results are shown in Fig. 4. For 
different stacked layer, the curves are almost parallel to each other as shown in Fig.4(a)~4(d). Only the curve for L=l 
is clearly differs from the others. In Fig. 4(e), the curves overlap with each other and can not be clearly separated. 
In Fig. 4(f), all curves start from the same point, but with different slope. The slop of other curves changes slightly 
except for the curve L=l. These results indicate that (i) the isotropic function for monolayer graphene is dramatically 
different from multi-layer graphene; (ii) from L=2 onwards, the isotropic function changes little with increasing the 
stacked layers. This is consistent with the results obtained from the information of the LDOS study that monolayer 
graphene is a special case in multi-layer graphene. Starting from bilayer, the electronic strcuture of graphene changes 
little when layers increase. 

V. CONCLUSION 

In summary, an analytical form of the real space Green's function (propagator) of multi-layer graphene is firstly 
constructed by using the effective-mass approximation. The Green's function demonstrates elegant spatial anisotropy 
with three-fold symmetry. Our numerical results show that the monolayer graphene is dramatically different from 
multi-layer graphene, and it is a special species in the family of graphite. In addition, we provides a new feasible 
method to identify the graphene layers by measuring the d 2 I/dV 2 in STM experiment, or the predicted features based 
on our simulated results can be used to verified STM measurements. Since multi-layer graphene is described by using 
the simple non-interactive tight-binding scheme, we are currently investigating how the Coulomb interaction plays a 
role in electronic structure of multi-layer graphene. 
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